/*
 * TotalEquilibriumSolver.cpp
 *
 *  Created on: 6 Dec 2011
 *      Author: allan
 */

#include "TotalEquilibriumSolver.h"

// GeoReact includes
#include <Core/CannonicalSystem.h>
#include <Utils/Assert.h>

TotalEquilibriumSolver::TotalEquilibriumSolver() 
{}

TotalEquilibriumSolver::TotalEquilibriumSolver(const CannonicalSystem& cannonical_system)
{
	Initialiaze(cannonical_system);
}

void TotalEquilibriumSolver::Initialiaze(const CannonicalSystem& cannonical_system)
{
	// Auxiliary variables
	const unsigned Nj = cannonical_system.GetNumPrimarySpecies();
	const unsigned Ne = cannonical_system.GetNumEquilibriumSpecies();
	const unsigned Nk = cannonical_system.GetNumKineticSpecies();
	const unsigned N  = Nj + Ne + Nk;
	
	// The indexes of the primary, equilibrium and kinetic species
	const vector<Index>& j_indexes = cannonical_system.GetIndexesPrimarySpecies();
	const vector<Index>& e_indexes = cannonical_system.GetIndexesEquilibriumSpecies();
	const vector<Index>& k_indexes = cannonical_system.GetIndexesKineticSpecies();
	
	// The canonical stoichiometric matrix of the equilibrium reactions
	const MatrixXd vej = cannonical_system.StoichiometricMatrixVej();
	
	// Assemble the equilibrium condition matrix (alpha)
	MatrixXd alpha = MatrixXd::Zero(Nj + Nk, N);
	
	for(unsigned j = 0; j < Nj; ++j)
	{
		alpha(j, j_indexes[j]) = 1.0;
		
		for(unsigned e = 0; e < Ne; ++e)
			alpha(j, e_indexes[e]) = vej(e, j);
	}
	
	for(unsigned k = 0; k < Nk; ++k)
		alpha(Nj + k, k_indexes[k]) = 1.0;
	
	equilibrium_conditions.SetConditionMatrix(alpha);
	
	equilibrium_solver.Initialiaze(cannonical_system.GetMultiphase(), 
		cannonical_system.GetEquilibriumReactions());
}

const ConvergenceDiagnostics
TotalEquilibriumSolver::Solve(const VectorXd& ujnk, double T, double P, VectorXd& n)
{
	equilibrium_conditions.SetConditionFunctions(ujnk);
	
	return equilibrium_solver.Solve(equilibrium_conditions, T, P, n);
}
